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Abstract. We give a goal-oriented a posteriori error estimator for the atomistic-continuum model- 
ing error in the quasicontinuum method, and we use this estimator to design an adaptive algorithm 
to compute a quantity of interest to a given tolerance by using a nearly minimal number of atom- 
istic degrees of freedom. We present computational results that demonstrate the effectiveness of 
our algorithm for a periodic array of dislocations described by a Frenkel-Kontorova type model. 



1. Introduction 

Multiscale methods offer the potentiaf to solve complex problems by utilizing a fine-scale model 
only in regions that require increased accuracy. However, it is usually not known a priori which 
region requires the increased accuracy, and an adaptive algorithm is needed to compute a given 
quantity of interest to a given tolerance with nearly optimal computational efficiency. 

The quasicontinuum (QC) method [3. 19I. Hol. [lit * s a multiscale computational method for crystals 
that retains sufficient accuracy of the atomistic model by utilizing a strain energy density obtained 
from the atomistic model by the Cauchy-Born rule in regions where the deformation is nearly 
uniform. The atomistic model is needed to accurately model the material behavior in regions of 
highly non-uniform deformations around defects such as dislocations and cracks. 

The approximation error within the quasicontinuum method can be decomposed into the mod- 
eling error which occurs when replacing the atomistic model by a continuum model, and the coars- 
ening error which arises from reducing the number of degrees of freedom within the continuum 
region. This paper purely focuses on the estimation of the modeling error. The optimal strategy 
to determine the mesh size in the continuum region will be studied in a forthcoming paper. 

The development of goal-oriented error estimators for mesh coarsening in the quasicontinuum 
method has been given in 0,0], and goal-oriented error estimators for atomistic-continuum mod- 
eling have recently been given in 0|. In all these works, the error is measured in terms of a user- 
definable quantity of interest instead of a global norm. Goal-oriented error estimation in general is 
based on duality techniques and has already been used for a variety of applications such as mesh 
refinement in finite element methods [l|, 0] and control of the modeling error in homogenization 0] . 

In 0], an a posteriori error estimator for the modeling error in the quasicontinuum method has 
been developed, analyzed, and applied to a one-dimensional Frenkel-Kontorova model of crystallo- 
graphic defects 0]. In this paper, we summarize this approach and adapt it to a different setting. 
Instead of clamped boundary conditions, we use periodic boundary conditions here which are phys- 
ically more realistic and allow for more succinct formulas. Moreover, an asymmetric quantity of 
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Figure 1: Frenkel-Kontorova model. The wells depict the misfit energy caused by the substrate. 



interest is used here rather than the symmetric one studied in [2| to give further insight into the 
behavior of the error estimator. 

The one-dimensional periodic Frenkel-Kontorova model chosen here allows for an easy study 
of the error estimator and keeps the formulas short, but still exhibits enough of the features of 
multidimensional problems for a realistic study of the error estimator. In addition to the nearest- 
neighbor harmonic interactions between the atoms in the classical Frenkel-Kontorova model, we add 
next-nearest-neighbor harmonic interactions to obtain a non-trivial quasicontinuum approximation. 

We describe the atomistic model and its quasicontinuum approximation in Section [21 and we 
formulate the error estimator in Section [3j We then develop in Section [3] an algorithm which employs 
the error estimator for an adaptive strategy, and we conclude by exhibiting and interpreting the 
results of our numerical experiments. 

2. Atomistic Model and Quasicontinuum Model 

As an application for the method of error estimation described in this paper, we study a periodic 
array of dislocations in a single crystal. We employ a Frenkel-Kontorova type model |5y to give a 
simplified one-dimensional description of these typically two-dimensional or three-dimensional phe- 
nomena. To model the elastic interactions, we consider 2M atoms in a periodic chain that interact 
by Hookean nearest-neighbor and next-nearest-neighbor springs. The dislocation is modeled by the 
interaction with a substrate which gives rise to a misfit energy, see Figure [TJ 

The vector y = (y^M+i, ■ ■ ■ , Um) £ M. 2M describes the positions of 2M atoms which generate the 
positions of an infinite chain of atoms by the relation 

Ui+2M = Vi + (2M + l)a for i = -oo, . . . , oo, (2.1) 

where ao denotes the lattice constant. The relation (|2.ip gives an average strain of (2M + 1)/2M 
due to a periodic array of dislocations that stretch the chain by one lattice constant every 2M 
atoms. 

The total energy 8 a for this atomistic system reads 



£ a (y) = £ a ' e (y) + £ a ' m (y) 



(2.2) 
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with elastic energy given by (recall (|2.1 
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i=-M+l i=-M+2 
M-l 

E (Vi+i ~Vi- ao) 2 + {y-M+i - VM + 2Ma f 
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i=-M+2 



+ {y-M+2 ~ VM + (2M - l)oo) 



and misfit energy given by 



£ a,m (y) 



M 

£ 

=-M+l 



Vi - ao 



a 2_ 



(2.4) 



Here ko, k% and ki denote the elastic constants. In the misfit energy, [^J denotes the largest 
integer smaller than or equal to x. We can obtain the following more symmetric form of the elastic 
energy (|2.3j) by realizing that the forces of constraint corresponding to the strain (|2.1p move the 
equilibrium spacing of the chain to 2 ^ M " 1 ao : 



+ \k 2 



M 



1 „2 



E [Vi+l -Vi~ ((» + 1) mod 2M) - i)o ) + 517",, 

i=-Atf+l 
M 



Li=-M+1 



2M+1 
2M 



(((i + 2) mod2M)-i)a + ^ag 



(2.5) 



where all atom indices are understood modulo 2M, and % mod 2M denotes i modulo 2M within 
the interval — M + 1, . . . , M. In the following, we neglect the constant terms since they do not have 
any effect when finding energy minimizers later. 

We consider a single vacancy between the atoms yo and y±. If we assume that the M leftmost 
atoms yi for — M + 1 < i < are situated in the interval ((i — |) oq, (i — ^) ao) and that the M 
rightmost atoms yi for 1 < i < M are situated in the interval ((i— 5) clq, (i + ^) ao), then the 
misfit energy can be expressed more simply as 



M 



£ a ' m (y) 



with 



-M+l 



(i — l)ao for i < 0, 
iao for z > 1. 



(2.6) 



To reduce the work in computing (|2.2p . we employ the quasicontinuum method @, [lO, 11] which 
has been successfully used for many applications. The quasicontinuum method consists of two steps: 
the passage to a continuum energy within the continuum region of the chain, and a subsequent 
coarsening in the continuum region to reduce the number of degrees of freedom. 
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In the first step, we replace the atomistic energy of all atoms from the continuum region by a 
continuum energy. To this end, we split the total energy into atom-wise contributions: 



M 
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The corresponding continuum energy can be derived following [2] to be 
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i=-M+l 



(2.7) 



(2.8 



(2.9) 



(2.10) 



(2.11) 



denotes the mixed atomistic-continuum energy. 

In the second step of the quasicontinuum approximation, the chain is coarsened in the contin- 
uum region by choosing representative atoms, more briefly called repatoms. The chain is then 
fully modeled in terms of the repatoms. The missing atoms are implicitly reconstructed by linear 
interpolation according to the Cauchy-Born hypothesis. The lengthy expression of the resulting 
quasicontinuum energy 



£ qc {y) 



(2.12) 



is not needed in this paper since we focus on the estimation of the modeling error. Hence we refer 
to [2] for the formula and its derivation. The error arising from coarsening will be studied in a 
forthcoming paper. 

For the subsequent argumentation, it is useful to rewrite the energies in matrix notation. We 
have 



£ a {y) = Uy - * a fD T E a D(y - a a ) + i(y - b a ) T K a (y - b a ) 



£ ac (y) = |(y - 8L a ) T D T E ac D(y - a a ) + |(y - h a ) T K a ' 



b a ), 



(2.13) 
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where the 2M x 2M matrices are given by 

D{ } i = — 1, = 1, 

(E a ) i4 = h + 2k 2 , {E a ) iA+1 = (E a ) i+1>i = k 2 , 

{E a %, = lh (6f + 5? +1 ) + \k 2 {8U + % + 5? +1 + Sf +2 ) + {\k x + 2k 2 ) {81 + 8f +1 ) , (2.14) 
(E a % l+1 = (E-% +l , = \k 2 (6? + 5f +2 ) , 
(K a ) ht = ko, 

with i = —M + 1, ...,M and all indices to be understood modulo 2M as before. The vectors 
a a e R 2M and b a G R 2M are defined as 

a*=[(-M + 1)3^00 (-M + 2)^ao ••• (M-l)^a M^a ] T , 

b a = [6-M+l &-A/+2 • • • &M-1 &Af] • 

We require that the elastic moduli satisfy k\ + 2k 2 > 2\k 2 \ to ensure that -E a is positive definite 
and that the misfit modulus k$ > to ensure that K a is positive definite. 

We are interested in finding energy minimizing configurations y a , y ac , and y qc of £ a , £ ac , and 
£ gc , respectively. The minimizers y a and y ac satisfy the linear equations 

M a y a = f a , 

(2.16) 

where 

M a := D T E a D + K a , f a := D T E a Da. a + K a b a , 

(2-17) 

M ac ._ D T E ac D + R a ^ f ac ._ D T E ac Dg a + R a h a 

We refer to (|2.16p as the primal problems. Note that the minimizers are uniquely determined due 
to the convexity of the energy. 

3. Error Estimation 

In the preceding section, we described how the quasicontinuum method gives an approximation 
y ac of the atomistic solution y a by passing from the fully atomistic model to a mixed atomistic- 
continuum formulation, and then we briefly mentioned how a further approximation, y qc , can be 
obtained by coarsening in the continuum region. 

Instead of measuring the error in some global norm, we measure the error of a quantity of interest 
denoted by Q(y) for some function Q : M. 2M — > R. We assume that Q is linear and thus has a 
representation 

Q(y) = q T y (3.1) 

for some vector q £ M. 2AI . We then have the splitting 

\Q(y a ) - Q(y qc )\ = \Q(y a - y ac ) + Q(y ac - y qc )\ < \Q(y a - y ac )\ + \Q(y ac - y qc )\ (3.2) 

of the total error into the modeling error, \Q(y a — y ac )|, and the coarsening error, |Q(y ac — y 9C )|, 
everything in terms of the quantity of interest. In this paper, we restrict ourselves to the estimation 
of the modeling error. The coarsening error will be analyzed in a forthcoming paper. 

An important tool for the error estimation in terms of a quantity of interest are the dual problems 
for the influence or generalized Green's functions g a and g ac given by 

M a g a = q, 

M ac g ac = ^ (3-3) 
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The matrices M a and M ac are symmetric since they stem from an energy, and we thus do not need 
to use their transpose for the dual problems. 

We denote the errors and the residuals, both for the deformation y a and the influence function 
g a , by 

e := y a - y ac , R(y) := M a (y a - y) = f° - M a y, 

e := g a - g ac , %) ■= M a (g a - g) = q - M a g. ( "' ' ' 

Then we have the basic identity for the error of the quantity of interest 

Q(y a ) - Q(y ac ) = q T e = g aT M a e = (g ac + e) T M a e 

= g acT R{y ac )+e T M a e. ( ' > " >) 

The quantities y ac and g ac are considered to be computable since the continuum degrees of 
freedom give local interactions, whereas y a and g a are not considered to be computable since they 
require a full atomistic computation. Thus the first term g acT R(y ac ) is easily computable, and the 
challenge is to estimate e T M a e. Let us note that in applications to mesh refinement for linear finite 
elements, the residual term vanishes due to Galerkin orthogonality, whereas in other applications 
it can be dominant over the second term. 

We utilize two error estimators derived in [2j and briefly summarized here. Our first error 
estimator is based on the generalized parallelogram identity 

e T M a e = l\\ae + <J~ l ef Ma - \\\ae - (r _1 e||^ (3.6) 

for all a ^ 0, where the M a -norm of some vector z is defined by ||z||jifo := (z T M a z) 1//2 . We define 
the computable bounds 

r?io W < Ike ± o-^ellAfa < rj± pp (3.7) 

by 

r?± p :=\\PD[a(y ac - a a )±a' 1 g ac ]\\ Ea , 

± __ \(y ac + e ± g ac ) T r ± \ ( 3 - 8 ) 
Vlow := \\y ac + d^hi- 

where 

P :=I -(E a )~ l E ac , 

(3.9) 

r ± := aR{y ac ) ±a~ 1 R(g ac ). 

Optimization of the bounds with respect to a and 6 leads in 0] to the following choice of the 
parameters: 



a : 



\PDg° 



i/ \\PD(y ac - aL a )\\ E a' 

q± r y g m y ~ r g iiy iim° 

yztTo-ac cfdcT ]\/fa-. r ac r ±T-. r ac ||rr ac ll^ 

1 & & ivi y 1 y n& iiM a 

Theorem 3.1. We have that 

\Q(y a )-Q(y ac )\<vi, (3.11) 

where the computable error estimator rj\ is defined as 

m := max {\g acT R(y ac ) + ±(r,+J 2 - \(^- pp ) 2 \ , \g acT R(y ac ) + I^) 2 - \{vT 0W ?\) ■ (3-12) 
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el 



We also developed the following weaker estimator in [2] using the Cauchy-Schwarz inequality 
in place of the parallelogram identity in ()3.6|) . We note that this estimator can be decomposed 
among the degrees of freedom and can thus be utilized in adaptive atomistic-continuum modeling 
decisions. 

Theorem 3.2. We have that 

M M 

\Q{y a )-Q{y ac )\<m< £ £ (3.13) 

i=-M+l i=-M+l 

where the computable global error estimator 7/2 and the computable local error estimators r/f* and 
r/l^, associated with atoms and elements, respectively, are defined as 

% := |g acT i?(y ac )| + \\PD(y ac - a a )\\ E a\\PDg ac \\ E a, 

\grR(y ac h\, i = -M + l,...,M, 

\ \(PD(y ac - a a )) .{{E a - E ac )D{y ac - a a )).| 
+5 \(PDg ac )i((E a - E ac )Dg ac ) i \ , i = -M + 1,...,M. 

4. Numerics 

Now we use the two a posteriori error estimators given in Section [3] to formulate an algorithm 
which adaptively decides between atomistic and continuum modeling. Finally, we present and dis- 
cuss the numerical results for the periodic array of dislocations described by the Frenkel-Kontorova 
model. 

The error estimator 771 gives a better estimate of the error because 7/2 involves additional inequal- 
ities. However, r/2 allows for an atom- wise decomposition, whereas r\\ does not. This is due to the 
fact that {f]f ow ) 2 in the definition of r/i is equal to the square of a sum of atom-wise components 
and not the sum of the square of these components. We can thus let r/i decide whether a given 
global tolerance r g i is already achieved or not and use the decomposition of r/2 to decide where the 
atomistic model is needed for a better approximation. In this way, we combine the better efficiency 
of r/i with the error localization of r/2. 

We start with a fully continuum model. We then switch to the atomistic model wherever the 
local error exceeds an atom-wise error tolerance r a t- While decreasing r a t, the algorithm adaptively 
tags larger and larger regions atomistic until the estimate for the goal-oriented error finally reaches 
T g i. The complete algorithm reads: 

(1) Choose T g i. Model all atoms as a continuum. Set T at <— T g i. 

(2) Solve primal problem (|2.16p for y ac and dual problem f)3.3|) for g ac . 

(3) Compute error estimator r/i from (|3.12j) . 

(4) If r/i < T g i, then stop. 

(5) Compute local error estimators r/^ and 77^ from (|3.14p . 
6 Set T at <- — • 

(7) Make all atoms i atomistic (5f = 1 and 5f = 0) for which 

r% ■■= + \ {vti-l + Vt) > rat- (4.1) 

(8) Go to (2). 

The factor t^v > 1 describes the rate at which the atom-wise tolerance T a t is decreased during the 
adaptive process. We found that Tdiv = 10 gives an efficient method for this problem. 
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iteration 


atomistic region 


Tat 


m 


1 


none 




1.000000e-10 


6.860546e-03 


2 


-26 . . . 


55 


1.000000e-ll 


1.238016e-07 


3 


-30 . . . 


60 


1.000000e-12 


2.600112e-08 


4 


-34 . . . 


66 


1.000000e-13 


3.922946e-09 


5 


-38 . . . 


73 


1.000000e-14 


4.104868e-10 


6 


-43 . . . 


80 


1.000000e-15 


4.105166e-ll 



Table 1: Convergence of the algorithm for r g i = 10 




-500 -400 -300 -200 -100 100 200 300 400 500 -500 -400 -300 -200 -100 100 200 300 400 500 



Figure 2: Decomposition of the error estimator vffl f° r iteration 1 (left, fully continuum model) and for 
iteration 6 (right, atomistic region —43 . . . 80). 



Now we come to the results for the Frenkel-Kontorova dislocation model for a periodic chain of 
1000 atoms, that is M = 500. The elastic constants are set to be /co = 1 and k\ = &2 = 2. To define 
the quantity of interest, we choose the average displacement of atoms 11 . . . 30. This leads to 

q = (qi)i=-M+l,...,M> U = 1 for 11 < i < 30, % = otherwise. (4.2) 

The global tolerance is chosen to be r g i = 10~ 10 . 

Table [T] shows how the successive adaptive determination of the atomistic-continuum modeling 
proceeds. After six iterations, the atom- wise tolerance is small enough so that 771 < r g i, that is the 
desired accuracy has been achieved. 

Figure [2] shows the decomposition rf^l of the error estimator 772 for the fully continuum model 
in iteration 1 of the adaption process. One can clearly read off from the graph that the error is 
large near the dislocation between atoms and 1 and near atoms 11 and 30, and that it decays 
exponentially away from these points. We note the slight nonsymmetry of the atomistic-continuum 
modeling due to using a goal function which averages over atoms 1 1 ... 30 to the right of the 
dislocation, but not to its left. The graph on the right shows the decomposition of r]2 in the final 
iteration 6 with an atomistic region given by indices —43 . . . 80. It exhibits the same nonsymmetry, 
but the error is considerably smaller with peaks at the boundary between the atomistic region and 
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atomistic 


IQ(y a -y ac )l 








V2 


rej 


$ion 


m 


\Q(y a -y ac )\ 


m 


\Q(y a -y ac )| 


none 


1.416421e-03 


6.860545e-03 


4.843577 


1.231314e-02 


8.693133 


-4 


.. 10 


1.863104e-03 


6.107510e-03 


3.278136 


1.049800e-02 


5.634680 


-9 


.. 20 


1.000572e-05 


3.358722e-04 


33.56803 


6.621488e-04 


66.17705 


-14 


. . 30 


1.430363e-04 


3.187552e-04 


2.228492 


5.140285e-04 


3.593694 


-19 


.. 40 


1.675490e-05 


2.626711e-05 


1.567727 


3.691344e-05 


2.203142 


-24 


.. 50 


7.361419e-07 


1.190138e-06 


1.616723 


1.693910e-06 


2.301065 


-29 


.. 60 


3.139276e-08 


5.157753e-08 


1.642975 


7.388556e-08 


2.353586 


-34 


.. 70 


1.146997e-09 


2.001550e-09 


1.745035 


2.934377e-09 


2.558312 



Table 2: Efficiency of the error estimators, ?7i/|Q(y a — y ac )| and r]2/\Q(y a — y oc )|- 



the continuum region. In both diagrams, the fluctuations come from the limited relative machine 
precision of about 10 -16 . 

Finally, Table [2] shows the efficiency of the error estimators 771 and 772 for different atomistic 
regions. |Q(y a — y ac )| gives the actual error which can be computed for this relatively small 
problem. In real applications, it is of course not available. One can clearly see that 7/1 gives a better 
estimate than 772, which numerically confirms our conjecture that 771 is a better estimator than 772. 
An unusually high value for the efficiency occurs when the atomistic-continuum boundary sweeps 
through the region where the quantity of interest is measured. After this, the efficiencies converge 
to decent values around 1.7 and 2.5 for 771 and 772, respectively. We note that for clamped boundary 
conditions and a symmetric quantity of interest, better efficiencies of 1.4 and 2, respectively, have 
been obtained [2|. 
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